Controlling the transverse instability of dark solitons 
and nucleation of vortices by a potential barrier 
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We study possibilities to suppress the transverse modulational instability (MI) of dark-soliton 
stripes in two-dimensional (2D) Bose-Einstein condensates (BECs) and self-defocusing bulk optical 
waveguides by means of quasi-lD structures. Adding an external repulsive barrier potential (which 
can be induced in BEC by a laser sheet, or by an embedded plate in optics), we demonstrate that it 
is possible to reduce the MI wavenumber band, and even render the dark-soliton stripe completely 
stable. Using this method, we demonstrate the control of the number of vortex pairs nucleated 
by each spatial period of the modulational perturbation. By means of the perturbation theory, we 
predict the number of the nucleated vortices per unit length. The analytical results are corroborated 
by the numerical computation of eigenmodes of small perturbations, as well as by direct simulations 
of the underlying Gross-Pitaevskii/nonlinear Schrodinger equation. 



I. INTRODUCTION 

Since the experimental creation of atomic Bose- 
Einstein condensates (BECs) [l|, a great deal of exper- 
imental and theoretical efforts has been invested into 
studies of nonlinear coherent matter waves. The intrin- 
sic nonlinear nature of BECs, originating from the inter- 
atomic interactions (accounted for by an effective mean 
field), together with the techniques available for the gen- 
eration of diverse initial configurations in the conden- 
sates, and the tunability of both external trapping po- 
tentials and intrinsic nonlincarity, have made it possible 
to study an extremely rich variety of nonlinear macro- 
scopic excitations [li, il]- In this context, it is relevant to 
stress that the effective dimensionality may be chosen as 
corresponding to the underlying three-dimensional (3D) 
geometry, or reduced to nearly 2D (pancake-shaped) and 
ID (cigar-shaped) settings, by applying a strong confine- 
ment in the "undesirable" direction(s) 0- On the other 
hand, the s-wave scattering length, which determines the 
nonlincarity strength in the framework of the mean-field 
description of the BEC, may be manipulated by dint of 
the magnetically- optically- 0] or confinement- Q 
induced Feshbach resonances. In combination with the 
use of properly designed trapping (magnetic and/or opti- 
cal) potentials, such as optical lattices 0, this technique 
paves the way towards a versatile and exceptionally ac- 
curate control over the nonlinear matter waves. Among 
so generated nonlinear modes, matter-wave solitons have 
been observed in a series of famous experiments. In par- 
ticular, nearly-lD bright and dark solitons were created 
in BECs with, respectively, attractive [sl-fioj and repul- 
sive pH - flTl interatomic interactions. Bright solitons of 
the gap type have also been created in the quasi-lD re- 
pulsive BEC [l^ loaded into an optical lattice. 



The description of the dynamics of matter waves is 
based on the Gross-Pitaevskii equation (GPE) which 
bears significant similarities to the nonlinear Schrodinger 
equation (NLSE) governing the transmission of light sig- 
nals in nonlinear optical media. Accordingly, the matter- 
wave solitons are counterparts of optical solitons, that 
have been studied in detail, in temporal, spatial |19| and 
spatio-temporal [1^ settings alike. 

In this work, we focus on dark solitons (see the re- 
views [U, for optical and matter- wave dark solitons, 
respectively), in connection to higher-dimensional geom- 
etry. In that case, an important issue is the stability of 
dark matter waves, which was first analyzed in the frame- 
work of the (2-|-l)-dimensional NLSE. In particular, the 
stability of the dark-soliton stripe (DSS) was studied in 
Ref. [23[ (see also Refs. [13, H^), where it was shown that 
the DSS is prone to the transverse modulational instabil- 
ity (MI) (alias "snaking instability") against transverse 
long-wavelength perturbations. Experimental [2^ 
and theoretical [2ll, [2^ studies of this instability in the 
context of nonlinear optics have revealed that it may lead 
to splitting of the DSS into a chain of vortices with alter- 
nating topological charges (vortex-anti- vortex pairs) . In 
particular, it was found that a quiescent ("black" )DSS is 
vulnerable to transverse "snaking" deformations [2ll [28| , 
causing the splitting into vortex pairs. Unstable moving 
("gray") DSSs do not split into vortices, but rather emit 
radiation in the form of sound waves. In the BEC con- 
text, the transverse MI and splitting of DSSs into vortex 
rings was first observed in the experiment with a two- 
component BEC composed of ^""Rb atoms in two different 
hyperfine states [2^. In that work, a DSS was created in 
one component, and the snaking instability caused it to 
decay into vortex rings (in a quasi-spherical geometry), 
in accordance with the theoretical predictions [s^l ■ 
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Several works proposed possibilities of suppressing the 
snaking instability of DSSs. In particular, in the context 
of nonlinear optics it is known [2ll | that this instab ility 



can be avoided in finite-size holding optical beams [31 



Simil arly , in the BEC context, it was shown [32[ (see also 
Ref . [33| ) that the snaking instability may be suppressed 
in sufficiently strong traps. On the other hand, it has 
been shown that 2D DSS can be stabilized by nonlocal 
nonlinearities [34|. Furthermore, 3D dark solitons [Hj 
may be stabilized in dipolar condensates, which are char- 
acterized by long-range dipole-dipole interactions (see a 
recent review of this topic in Ref. [1^). In that case, the 
suppression of the DSS instability is provided by the non- 
local character of the respective mean-field model (see 
also Ref. [s^l for similar results in the context of op- 
tics). A more complex dark-soliton configuration, which 
is not subject to the MI, was reported in Refs. [3^ 
in the context of the two-component BEC: it features a 
"cross" formed by the intersection of two rectilinear do- 
main walls, with the wave functions of the same species 
filling each pair of opposite quadrants, with a tt phase dif- 
ference between them. In this way, a quasi-dark-soliton 
configuration is formed, which is stable for long times 
(even in the rotating trap) in a large parametric region. 

In the present work, we propose and analyze an experi- 
mentally relevant scheme for suppressing and controlling 
the transverse instability of DSSs in BECs and optics, 
based on the use of a repulsive quasi-lD potential barrier, 
which can be generated by a blue-detuned laser beam in 
BECs (see, e.g., Ref. (i^), or by a slab with a lower value 
of the refractive index embedded into the self-defocusing 
bulk waveguide. This barrier, which repels the atoms in 
the condensate or the light in the waveguide, creates a 
channel where a DSS can be naturally trapped. Ana- 
lyzing the transverse MI of the DSS in this setting, we 
show that the channel can suppress the snaking insta- 
bility. In fact, results produced by our linear-stability 
analysis and corroborated by direct simulations indicate 
that the wave number of the most unstable perturbation 
mode is shifted due to the presence of the potential bar- 
rier. The control of the most unstable mode is useful not 
only for the complete stabilization of the DSS, but also 
for controlling the number of vortex pairs generated by 
the snaking instability. This control mechanism can be 
adjusted by appropriately tuning the height and trans- 
verse width of the potential barrier. 

The paper is organized as follows. In Sec. |IT1 we intro- 
duce the model and produce an approximate solution for 
the DSS, in the presence of the potential barrier. Using 
the variational approximation, we obtain a solution for 
the DSS and analytically predict the critical wavenum- 
ber for the MI modes, by means of the linear-stability 
analysis (details for the latter are described in the ap- 
pendix). In Sec. Illli we numerically study the trans- 
verse MI of the DSS, varying the height and width of the 
stabilizing potential barrier. We find that, in the uni- 
form space (i.e., in the absence of an external harmonic 
trapping potential) and with a sufficiently broad repul- 



sive potential barrier, there exists a critical value of its 
strength (height), above which DSSs get completely sta- 
bilized against the MI. In the case when the instability 
is not completely suppressed, we explore the possibility 
of controlling the density of the nucleated vortex pairs 
for an infinitely long DSS, by tuning the most unstable 
mode, via the parameters of the barrier. We then con- 
sider the more realistic model of a harmonically confined 
BEC (which may be relevant to optics too, representing 
a rod-shaped bulk waveguide), and perform numerical 
simulations which demonstrate the degree of the control 
over the number of the nucleated vortex pairs. In Sec. lIVl 
we summarize the findings and point out directions for 
future work. 



II. THE MODEL AND ITS ANALYTICAL 
CONSIDERATION 



The fundamental equation and dark-soliton 
solution 



We start by considering the following (2-1-1)- 
dimensional GPE/NLSE in the usual scaled form, with 
the repulsive nonlinearity: 

In the BEC context, this equation governs the evolu- 
tion of the macroscopic wave function u{x,yA) of the 
pancake-shaped condensate in the (x, y) plane [l,!!]. The 
external potential is assumed to have the following form: 



Vix,y) ^VHTix,y) + Vi^six), 



(2) 



that includes an external harmonic trap with frequencies 



VHT{x,y) 



(3) 



and the potential barrier (corresponding to the far- 
detuned laser sheet illuminating the BEC) with the 
Gaussian profile: 



Fls(x) = Aexp[-xV(2a2)] . 



(4) 



Here A and a measure the height and width of the bar- 
rier, with A > Q {A < 0) corresponding to the blue- (red-) 
detuned laser beam, which repels (attracts) atoms in the 
condensate (note that yl < corresponds to a potential 
trough of depth \A\, rather than a barrier). 

In terms of nonlinear optics, Eq. ^ governs the evo- 
lution of the local amplitude of the electromagnetic wave 
in the self-defocusing bulk waveguide, with t replaced by 
the propagation distance z, while x,y are the transverse 
coordinates, and V{x, y) describes a local modulation of 
the refractive index. The barrier potential ^ then cor- 
responds to a slab with a lower value of the refractive 
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FIG. 1: (Color online) Integral (|12p as a function of the 
rescaled width a = /3a. 



index embedded into the waveguide, while potential ^ 
approximates for the global waveguiding structure. 

In the absence of the external potential, V{x,y) = 0, 
Eq. ([!]) admits an exact analytical solution for the DSS. 
Assuming that its nodal plane is oriented along the x- 
dircction, to be aligned with the potential barrier, when 
it is switched on, the DSS solution is 



u{x, y, t) = y/JI{Btanh[y/JIB{x — vt)] + iA} 



(5) 



Here, y^exp(— i^f) is the stationary background with 
density fj, (which is actually equal to the normalized 
chemical potential) supporting the dark soliton. Further, 
and v = ^/Jj^A represent, respectively, the soliton's 
depth (or inverse width) and velocity, with A and B sub- 
ject to the constraint A^+B^ = 1. The DSS with A^O 
{A 7^ 0) is stationary (moving), and is usually called 
a "black" ("grey") dark soliton. Below, we first employ 
the variational approximation to find an approximate sta- 
tionary DSS solution to Eq. ([T]) in the presence of the po- 
tential barrier (j?]), and then perform the linear-stability 
analysis to study the transverse MI. 



B. Dark-soliton stripe steady state in the presence 
of the potential barrier. 

In the absence of the harmonic trapping potential ([3|) , 
but in the presence of barrier (j4]) , the profile of a station- 
ary DSS solution to Eq. ([1]) is sought for as 



uix,y,t) = Uix)e-'^', 



(6) 



where ^ is the chemical potential (or — /i is the propaga- 
tion constant, in terms of the optical model), and U{x) 
is a real function. Substituting this into Eq. ([IJ leads to 
an ordinary differential equation. 



(7) 



As Eq. ([7]) describes both the dark soliton and the uni- 
form background, one should subtract the contribution 
of the latter into the respective Lagrangian, which will 
be used as the basis of the variational approximation. To 
this end, we follow Rcf. [4l| (see also Ref. [l^), and write 



the Lagrangian as follows. 



{U'^fi) + {U'-^i'), (8) 



where fi is the amplitude of the uniform background, 
as defined above. Then, to approximate the station- 
ary dark-soliton solution of Eq. ([7]) we use the following 
ansatz [cf. Eq. dS])] 



U{x;A,(t) = y/]Itanh{(3{A,a)x), 



(9) 



where the inverse width of the soliton, /3, is a varia- 
tional parameter, while the density fi is the fixed back- 
ground chemical potential. Inserting ansatz (|9]) into 
the Lagrangian ([5]) and performing the integration over 
—oo<x< +00, wc arrive at the effective (averaged) 
Lagrangian: 



+ CXD 



-^'/(2'^')sech2(/3a;)da;. 

(10) 

Further, we use the Euler-Lagrange equation, dLcs/d/S — 
0, to derive the following implicit equation for the inverse 
soliton's width (3: 



=M-3A/3V(/3,a), 



(11) 



where 



+ 00 



J(/3,ct)= / xe-'''/(^''')sech2(/3a;)tanh(/3a;)da;. (12) 



Using the rescaling, x — >■ /3x, we conclude that J(/3, <t) = 
J(l, ct)//3^, where = /3(t, hence one can rewrite Eq. (fTTj) 
as 



=H- 3AJ(l,/3cr). 



(13) 



Apparently, Eq. (|13p . with J being a function of the sin- 
gle argument, which is depicted in Fig. [T] versus rescaled 
width CT, is simpler than Eq. (jlip . 

In Fig. [5] we validate the results of the variational ap- 
proximation by comparing the width extracted from the 
numerically found steady state, and from the variational 
equation ([T^ [thick (blue) and thin (red) lines, respec- 
tively], for different values of the barrier's parameters. As 
seen in the figure, for ^ > (i.e., for the repulsive poten- 
tial barrier) , the variational approximation is very accu- 
rate, with the steady-state solution closely resembling the 
tanh profile of ansatz ^ [see panel (b) in the figure]. On 
the other hand, it is seen in panel (c) that for yl < (the 
attractive potential trough, instead of the barrier) the 
approximation deteriorates, especially for large widths 
of the trough (cr > 0.6, which exceeds the width of the 
DSS). This observation is explained by the fact that the 
steady-state profile develops a localized hump induced 
by the attractive trough (see also Ref. H^), which is not 
captured by our ansatz. 



CO. 
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FIG. 2: (Color online) (a) Inverse width /3 of the dark soliton 
as a function of width a of the potential barrier, with height A 
ranging from A = —1 (top curves) to A = -1-1 (bottom curves) 
in steps of 0.25. The thick (blue) curves depict the width 
of the numerical solution, obtained by fitting its profile to 
tanh(/3a;). The thin (red) lines represent the results predicted 
by the variational approximation, see Eq. (|13(l . Panels (b) 
and (c) depict, respectively, two solutions for {A, a) = (1,1) 
and {A, a) = (—1,1.5), the latter case corresponding to the 
potential trough, rather than a barrier. The thick (blue) line 
depicts the numerically found steady state, while the thin 
(red) line represents the fitted tanh(/3a;) profile. The trans- 
verse profile of the potential barrier is depicted by the dashed 
(green) line. The chemical potential is fixed to be ^ = 1. 

C. The transverse modulational instability of the 
dark-soliton stripe in the presence of the potential 
barrier. 

We now consider the solution of Eq. ([T|) in the form of 
the DSS in the presence of the potential barrier. In this 
case, the steady state DSS may be approximated by 

uix, 2/, t) = uo{x, y) e-*^* = ^UnhilSx) e"*^*, (14) 

where /3 is determined by Eq. (fT3|) . To analyze the MI 
of this DSS, we follow the approach of Ref. [231 and con- 
sider small transverse perturbations, see details in the 
appendix. The result of the analysis presented in the ap- 
pendix is that the critical wavenumber which defines the 
width of the instability band for solution ^ is 

fccr = /3, (15) 

where (3 is determined by Eq. ([13]). That is, the DSS 
p4[) is unstable against transverse perturbations with 
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FIG. 3: (Color online) The dynamics induced by the trans- 
verse modulational instability of a dark-soliton stripe in the 
absence of the external harmonic trap. Shown are different 
density snapshots at indicated moments of time. The initial 
condition corresponds to a stationary stripe with an added 
random perturbation of relative size 10~^. The top panel: 
the case of the attractive potential trough, with A = —1 and 
a = 1.5. The middle panel: no potential barrier or trough 
{A = Q). The bottom panel: the case of the repulsive poten- 
tial barrier, with A = 1 and a — 1.5. Note the faster destabi- 
lization and nucleation of a larger number of vortices (see dark 
circular spots for later times) in the case of the attractive po- 
tential trough (the top panel) , in comparison to the case when 
the trough is absent (the middle panel), and the complete sta- 
bilization of the dark-soliton stripe by the repulsive potential 
barrier in the bottom panel. The simulations in this figure 
were carried out in the domain of {x, y) G [—15, 15] x [—15, 15]. 



wavenumbers k < kcx =^ /3. When this instability sets 
in, the DSS undergoes a snake-like transverse deforma- 
tion, which eventually breaks it into vortex-antivortex 
pairs. This outcome of the evolution can be observed in 
Fig. El 

From Eq. (1151) we can derive the critical number of 
vortices in the chain produced by the splitting of the 
DSS dm): 

N„ = (3L/7:, (16) 

where L is the length of the DSS. Since f3 is implicitly 
determined by height A and width a of the potential 
barrier (|4]), as per Eq. (fT3)) . it is possible to control the 
critical wavenumber by appropriately adjusting the bar- 
rier's parameters. This, in turn, offers a control of the 
upper bound of the vortex pairs that are nucleated by 
the transverse instability of the DSS. In the next section 
we show that it is indeed possible to control the number 
of the nucleated vortex pairs, which we compare with the 
analytical prediction. 
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FIG. 4: (Color online) The growth rate of the transverse mod- 
ulational instability of dark-soliton stripes at different values 
of the height and width of the potential barrier, A and a. The 
panels correspond to the increasing width: (a) a — 0.5, (b) 
(7 = 1, and (c) a — 1.5, as indicated by the dashed lines in 
Fig. [5] Each panel displays the real part of the eigenvalue, 
for the potential trough/barrier with heights A = —1 [(blue) 
circles], A — [(green) squares], and A = 1 [(red) triangles]. 
The computations corresponding to this figure were carried 
out in the domain of {x,y) G [—50,50] x [—50,50]. 



III. NUMERICAL RESULTS 

A. Transverse modulational instability of the 
dark-soliton stripe in the absence of the harmonic 
trap 

The main subject of this work is the control of the 
transverse MI of the DSS by means of the potential bar- 
rier. Typical examples of the snaking instability arc de- 
picted in Fig. |3l in the absence of the external harmonic 
trap [oJx,y = in Eq. ^]. The figure depicts the evolu- 
tion of the density \u{x,y,t)\'^ at different times, start- 
ing from a stationary DSS with a small initial random 
perturbation added to it. The middle panel in the fig- 
ure depicts the case of a quiescent (black) DSS (/? = 1), 
with no potential barrier A ~ 0, where the instability 
manifests itself with wavenumber k « 27r/10 = 0.628 
(wavelength ~ 10), which is smaller than the predicted 
critical wavenumber, /ccr = /3 = 1. In fact, this unsta- 
ble wavenumber is the most unstable one, fcmax, for this 
configuration. 

In order to quantify the instability and to com- 
pare it to the analytical results of the previous sec- 
tion, we compute the stability spectra of the stationary 
solutions uo{x,y) through the standard Bogoliubov-de 
Gennes (BdG) analysis. This analysis involves numer- 
ically solving the linear eigenvalue BdG problem, which 
stems from the linearization of the GPE ([T]) around the 
steady state solution uq by using ansatz u{x,y,t) — 
{uo{x,y) + [a{x,y)e^* + b*{x,y)e^**]}e~^''K The solu- 
tion of this BdG eigenvalue problem yields the eigenfunc- 
tions {a{x,y),b{x,y)} and the eigenvalues A. Note that, 
due to the Hamiltonian nature of the system, if A is an 
eigenvalue of the BdG spectrum, so are also —A, A* and 
— A*. Notice that the linear stability condition amounts 
to Rc(A) = 0, i.e., all eigenvalues must be imaginary. In 
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FIG. 5: (Color online) Regions of the transverse modula- 
tional instability of the dark-soliton stripe for different pa- 
rameters of the potential barrier. Depicted is the real part of 
the eigenvalue (the lighter color represents the stronger insta- 
bility, while black corresponds to the stability) as a function 
of wavenumber k and the barrier's height A. Different panels 
correspond to increasing width of the barrier: (a) a — 0.5, (b) 
a = 1, and (c) a = 1.5. The black and white solid lines de- 
pict the prediction of the critical wavenumber as per Eq. (|15|l . 
Vertical dashed lines represent the cuts depicted in Fig. H) 



fact, in connection to the stability analysis described in 
the previous section, the BdG eigenvalue A corresponds 
to eigenfrequency through relation f2 = iX. Therefore, 
an unstable eigenvalue corresponding to Re(A) > cor- 
responds, in turn, to an eigenfrequency with Im(f2) > 0. 

The MI spectrum for the DSS in the case of no poten- 
tial barrier (A = 0) is depicted by the (green) squares 
in all panels in Fig. |4l The most unstable wavenum- 
ber, fcijiax, corresponds to the location of the maxima of 
Re(A) > 0. Note that, as predicted by the linear stabil- 
ity analysis, all modes with k > kcr = 13 = 1 are stable. 
It is also interesting to note that for strong enough at- 
tractive (A < 0) potential barriers the whole spectrum is 
shifted to the right, so that small wavenumbers become 
stable [see (blue) circles in Fig.H]. Furthermore, it is im- 
portant to note that the results presented in Fig. 2] were 
obtained in a large box, {x,y) G [—50,50] x [—50,50], 
in order to safely capture small wavenumbers. In con- 
trast, the simulations depicted in Fig. [3] were performed 
in a smaller domain box, {x,y) G [—15,15] x [—15,15] 
so that the case for A = 1 (the bottom row in the fig- 
ure) is rendered stable since the corresponding unstable 
wavenumbers do not fit into the integration box. This 
effect, discussed further below, becomes important when 
considering trapped condensates that inherently possess 
a finite size. 

We now turn to the effect of the potential barrier on 
the stability (and dynamics) of the DSS. The top panel 
in Fig. [3] depicts the snaking instability corresponding 
to the attractive potential trough with strength (depth) 
A ~ —1. As seen in the figure, the attractive trough 
naturally enhances the instability, in comparison to the 
free-standing DSS, in two ways: (i) the MI sets in earlier, 
and (ii) the most unstable wavenumber is right-shifted 
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(fcmax ~ 27r/6 = 1.0472). The corresponding instability 
spectrum for this case is depicted by the (blue) circles in 
Fig-EJc). In lact, comparing the (green) squares and the 
(blue) dots in this panel, it is evident that the net elFect 
of the attractive potential trough is the right-shift of the 
instability spectrum. Note that the spectrum is located 
to the left of the threshold, fccr = /? ~ 1.78, which is 
predicted by the linear-stability analysis. 

In contrast to the destabilizing effects of the attrac- 
tive potential, the repulsive barrier stabilizes the DSS, 
as might be expected and is shown in the bottom panel 
of Fig. [3] The spectrum for this case is depicted by the 
(red) triangles in Fig. |4|c). Note that, although some 
perturbation modes are still unstable in this case, their 
wavelengths are too large to fit the integration domain 
used in the simulations shown in the bottom panel of 
Fig. [3] and, thus, the respective sufficiently strongly re- 
pulsive barrier renders the DSS effectively stable (in a 
still larger domain, we do observe the development of 
the weak transverse MI in this case, which is not shown 
here). The same stabilizing effect is also observed for 
other values of the widths of the barrier, as it can be 
clearly seen in all spectra depicted in Fig. This figure 
depicts, for three values of width ct, the shift of the MI 
spectrum for the attractive [(blue) circles], zero [(green) 
squares] and repulsive [(red) triangles] channel potentials. 
The figure suggests that, as mentioned above, for a suffi- 
ciently strong repulsive potential barrier, only very small 
wavcnumbers may be unstable. Therefore, if the longi- 
tudinal size of the condensate (set by the external trap) 
or optical waveguide, into which the DSS is embedded, is 
not too large, the instability might be fully suppressed, 
as the unstable wavelengths would be too large to comply 
with the longitudinal size (see also Ref. [s^ and Scc. lIII Bl 
below). 

In fact, as we now demonstrate, a sufficiently strong 
repulsive barrier may completely suppress all instabili- 
ties, including those with very small wavcnumbers. In 
Fig. [51 we depict the instability spectra for three differ- 
ent widths of the barrier widths, as a function of the sign 
and strength of the channel potential. The light-colored 
areas in the panels correspond to unstable modes, with 
the shading scale corresponding to the associated insta- 
bility growth rate for each mode. As shown in the figure, 
for sufficiently wide potential channels {a > 1.0, see the 
two right panels), there is a positive value of the bar- 
rier's height above which all wavcnumbers arc stable (the 
zero real part of the eigenvalues corresponds to the black 
shade in the panels). Therefore, it is possible to com- 
pletely suppress the snaking instability and render the 
DSS stable with the appropriate choice of the potential 
barrier. 

For example, for the domain used in Figs. |4][5] and a 
potential barrier width of cr = 1 (cr = 1.5) and height 
larger or equal to A = 1.65 {A = 1.35) the DSS is com- 
pletely stable. We also plot in Fig. [5] (see solid black 
and white lines) the critical wavenumber (|15|) obtained 
in the analytical form in Sec. [ill This approximate result 



(valid for A = 0) captures the qualitative behavior of the 
critical wavenumber, but fails to give its precise location. 

Another important feature of the control over the lo- 
cation of the most unstable mode fcmax is that it allows 
one to precisely manipulate the number of nucleated vor- 
tices per unit length which emerge from the snaking in- 
stability. The instability nucleates one vortex pair per 
snaking wavelength [s^. This feature is clearly visible 
at later times in Fig. [3l For example, the destabilizing 
(attractive) channel potential, with {a, A) = (1.5,-1), 
enhances the outcome by producing iV„ = 10 vortices in 
the domain of length L = 30, namely the vortex density 
Pv = 10/30 — 0.33... . On the other hand, the free- 
standing DSS {A = 0) is responsible for the nucleation of 
Pv = 6/30 = 0.2 vortices per unit length. Thus, one can 
manipulate the density of the nucleated vortices, select- 
ing it from zero to a maximum value which is inversely 
proportional to the healing length of the condensate, in 
the case of BEC (recall that two single-charged vortices 
cannot coexist at distances smaller than this length). 



B. The modulational instability of dark-soliton 
stripes in the harmonic trap 

We now consider the combined effects of the external 
harmonic trap and the channel potential. As mentioned 
above, this combination can produce an effective sup- 
pression of the snaking instability if the trap's strength 
is such that the corresponding Thomas-Fermi (TF) ra- 
dius of the BEC, i?TF = -s/S/I/a; (for uj^ = ujy = lu), 
is smaller than the smallest wavelength of unstable per- 
turbations. In Fig. [5] we display the development of the 
DSS' snaking instability and the concomitant vortex nu- 
cleation for the BEC confined in the harmonic trap of 
strength w = 0.05. The three panels correspond, from 
top to bottom, to channel potentials with width cr = 1.5 
and strength A ~ —0.8, 0, and 0.4. It is difficult to pre- 
cisely measure the number of nucleated vortices because 
some of them are nucleated at the rims of the configura- 
tion and are almost invisible. Therefore, we herein focus 
on measuring the number of nucleated vortices, Ny, in- 
side the circle of the TF radius. It is clear that, for 
A = —0.8, at least Ny = 16 vortices are nucleated, while 
Ny = 10 for A ~ (without the channel potential), and 
only A^^ = 5 vortices emerge for A = 0.4. 

Figure [7] shows the number of nucleated vortices for 
two different values of the trapping frequency. The curve 
corresponding to lu = 0.05 has been scaled by factor 1/2 
for a better comparison with the one corresponding to 
ll) = 0.1, since the former has the TF radius which is 
twice as large. In order to estimate analytically the num- 
ber of nucleated vortices inside the BEC cloud, we first 
establish an upper bound by considering the number of 
nucleated vortices in a DSS of length L = 2i?TF with 
the threshold wavenumber k = kcr obtained in Sec. HH 
Clearly, this estimation is merely an upper bound for the 
number of nucleated vortices because of the following rea- 
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FIG. 6; Dynamics of the transverse modulational instability 
of a dark-soliton stripe in the BEG loaded into the harmonic 
trap of strength to = 0.05. Shown are different snapshots of 
the density at the times indicated in the figure. The initial 
condition corresponds to a stationary stripe to which a ran- 
dom perturbation of relative size 10~^ was added. The top 
panel shows the case with the attractive potential trough of 
depth A — —0.8 and width a = 1.5. The middle panel: the 
situation in the absence of the channel potential {A = 0). The 
bottom panel: the case of the repulsive potential barrier of 
height A = 0.4 and width a = 1.5. Note that the number of 
nucleated vortices can be controUably varied from about 16 
for A = —0.8 (in the top panel) to none for j4 > 1 (not shown 
here) . 



sons: (i) since k^ax is not readily available from the the- 
ory, we are using fccr (^max < ^cr), which corresponds 
to the instability threshold and not to the wavenumbcr 
with the maximum growth rate; (ii) at the periphery of 
the TF cloud, it is difficult to pick the number of vor- 
tices through the direct count of the number of empty 
cores in the density profiles, and (iii) as it can be noticed 
from Fig. [5l the analytical prediction for fccr is always an 
overestimation of the true wavenumber threshold. Thus, 
an upper bound for the number of vortices can be safely 
given by 



(17) 



In connection with item (iii) above, it is worth mention- 
ing that it is in principle possible to measure the number 
of vortices by analyzing the vorticity (curl of the super- 
fluid velocity) of the condensate. This method reveals 
(results not shown here) that, indeed, some vortices ere- 
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FIG. 7: (Golor online) The number of nucleated vortices [N^) 
resulting from the snaking instability of the dark-soliton stripe 
in the BEG loaded into the harmonic trap, versus the strength 
[A) of the channel potential. The number of vortices is mea- 
sured by direct count of the number of empty cores present in 
the density profile after the snaking instability is fully devel- 
oped. The solid line depicts half of the theoretical prediction 
ofEq. (|17p (see text). Gircles: the trap's strength corresponds 
to oj = 0.05 and depicted is the half number of the vortices 
{Nv/2). Squares: the trap's strength corresponds to cj = 0.1. 
The width of potential channel is a = 1.5 for both data sets. 



ated at the periphery of the cloud are missed by the direct 
empty-core counting method, based on the inspection of 
the density profiles. Nonetheless, we have observed that 
these peripheral vortices are quickly pushed out and al- 
ways stay at the periphery. Therefore, such boundary 
vortices do not induce a considerable change in the long- 
term BEC dynamics. It is also worth pointing out that 
counting vortices by the number of empty cores in the 
density profile is tantamount to what is often used in 
current BEC experiments. Of course, in order to im- 
age the vortices in real experiments, the BEC cloud has 
to be left to expand and this could potentially modify 
the density distribution and "push" vortices inwards or 
outwards the BEC cloud. However, since the expansion 
time is usually short, it is likely that any dynamics dur- 
ing expansion will be rather slow and, thus, no significant 
change in the measurable number of vortices should be 
expected. 

As concluded above, the observable number of vortices 
nucleated by the DSS, N^,, will be smaller than the up- 
per bound A^cr defined in Eq. (|T7| , hence a relation of the 
type Ny = aNcr, with a < 1, should be more appropri- 
ate. For the specific value of a = 1.5 used in Fig. [3 we 
have noticed that the most unstable growth rate fcmax is 
very close to half oi the critical threshold fccr for values of 
the laser intensity A < 1. Therefore, for this regime, it is 
sensible to choose a = 0.5 as a relevant approximation. 
In fact, as it can be noticed from Fig. [71 sec solid black 
line, we have found that the choice of a = 0.5 gives a very 
good match to the actual number of observed vortices in 
the numerical experiments for a wide range of values of 
the potential barrier intensity A. It is important, how- 
ever, to stress that the above estimate, relying on the 
assumption that fcmax ~ 0.5fccr, is only valid for a wide 
potential barrier with cr = 1.5. For other regimes, we 
can only give an upper bound for the number of nucle- 
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ated vortices as per Eq. (|T7)) . Furthermore, it is relevant 
to mention that the prefactor a also depends (results not 
shown here), for other fixed values of the laser sheet's 
width cr, on the laser strength A. 

IV. CONCLUSIONS 

In this work, we have proposed and analyzed a tech- 
nique to control and even eliminate the transverse modu- 
lational instability, alias the snaking instability, of dark- 
soliton stripes in Bose-Einstein condensates (BECs) and 
bulk optical waveguides with the self-defocusing nonlin- 
earity. This technique relies on the use of a quasi- ID po- 
tential barrier, that may be induced by a blue-detuned 
laser sheet in the BEC, or by a slab embedded into the 
optical waveguide. We have also considered the opposite 
quasi-lD potential, in the form of an attractive poten- 
tial trough (that may be induced by a red-detuning laser 
beam) , which enhances the snaking instability by making 
the instability-onset time shorter and, more importantly, 
the unstable wavenumbers larger. Therefore, the attrac- 
tive trough produces more instability-induced oscillations 
per unit length, eventually generating a larger density of 
vortex- antivortex pairs. Most importantly, the repulsive 
barrier is able to suppress the snaking-instability band, 
and even completely stabilize the dark-soliton stripe. 

We have also developed an analytical approximation 
for estimating the largest unstable mode (in the absence 
of the external harmonic trap). The analytical results 
show good agreement with the numerical results obtained 
from the Gross-Pitaevskii equation. We also considered 
the model including the isotropic harmonic trap. In this 
case, we monitored the total number of vortices nucleated 
by the modulational instability inside the BEC cloud. 
Assuming that the size of the cloud is determined by 
the Thomas-Fermi radius, we were able to provide an 
upper bound for the estimate of the nucleated number of 
vortices as a function of the harmonic-trap's strength. 

The stabilization mechanism proposed in this work 
may be naturally extended in diverse directions. First, 
it is interesting to consider the stabilization of quasi-lD 
solitons by the attractive potential trough in the model 
with the self-attraction, and the situation with either sign 
of the nonlinearity, if the channel potentials are replaced 
by a sufficiently strong periodic optical lattice. Another 
challenging possibility is to consider the stabilization of 
quasi-lD solitons by the cigar-shaped potential and its 
periodic optical-lattice counterpart in the 3D geometry. 
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Appendix: The linear stability analysis for the 
dark-soliton stripe 

In order to analyze the linear stability of the dark- 
soliton stripe characterized by Eq. ([H]) in the presence 
of the potential barrier, we closely follow the approach of 
Ref. [231. We consider a perturbation of the DSS in the 
form of 

6u^ui{x)e^, 5u* — U2{x)e^ , (18) 

with A = ~iflt+iky, where * denotes the complex conju- 
gation, Q is an eigenvalue in the spectrum of Eq. ([T]) lin- 
earized around solution ([T4|) . (w 1,142) is the correspond- 
ing eigenfunction, and k is the corresponding wavenum- 
ber along the y direction. 

Substituting the standard ansatz for the perturbation 
around the soliton solution (|T4)) . 

u{x,y,t) = {uQ + 5u + Su*)e~''''' (19) 

into Eq. ([I} , and noting the linear independence of eA and 
, leads one to the following linear-stability equations: 

^UqU2 — hk^ui + riui = 0, 
< ^ (20) 

^|^ + (M-^-2KI>; 

— U^U\ — ^k^U2 — ^2^2 = 0, 

where the original second equation has been replaced by 
its complex-conjugate counterpart. Equation (|20p can be 
written in the matrix form, 

nMip-^k'^Iip + Lip = 0, (21) 

where 

r = ir^_/'2KP-M + ^ u'o \ 

~ 2 dx^ \ uf 2\uo\^-fi + V J ■ 
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Rotating the cigcnfunction Lp, 



1 i 



yields another representation of Eq. (HI]): 



1 i 

-1 i 



with (p : 



h 



Using the matrix 1/2 



1 1 



muhiply Eq. ([23]) yields 



VlMnp - ^k'^If + Lip = 0, 



(22) 



^ = 
(23) 
to left 



(24) 



with Mj = 



i 


where A = 2\uo\^ + - ^l + V and A± = {u^'^±ul)/2. 
By using Eq. (fH)) . L becomes 



A iA_ 
iA_ A 

2 _L „,2N 



2 



ii+ 
L2+ 



(25) 



where we define 

Li+ = 3Aisech^(/3a;) - 2/i - L; 



2+ 



/i sech^(/?a-) — K 
(26) 

As follows from Eqs. (|2T]) and ([24]) . the main piece of 
the information concerning the spectrum of eigenvalues 
CI comes from properties of the operator L. Obviously, 
it is a Hermitian operator, hence its spectrum is purely 
real. The respective eigenvalue problem can be written 
as 



(27) 



We should note that p corresponds to k in Eq. ([M]) with 
f2 = 0. Our purpose is to find the critical value kd- of 
wavenumber k for the instability of solution (|14p . To 
proceed, we need to look for a solution to ([?7)) as 



-Px 
-fix 



92{x) 



(28) 



Substituting Eq. ([28]) into Eq. ([27]) leads to the following 
equations for gi{x) and 32(2;): 



dgi ^ 1 d^9i 
dx 2 dx'^ 

+ - 2/_t + 3//sech^(/3a;) - F ) .gi = 0, 
da; 2 dx'^ 

-ip^ + ^/3^ + Msech2(/3x) - .92 = 0. 

(30) 

Considering the asymptotic behavior of Eq. ([5(1]) in the 
limit \x\ — >■ 00, noting that V — > and sech^(/3a;) — >■ 0, 
and taking Eq. (|^^ into regard, we obtain 



^P' + ^/3' - 2/i ) 51 = 0, 
+ ^/3^ ) .92 0, 



(31) 



where = gid ± oo|) and 52 = 52(| i oo|) are constants. 
Equation (|5T]) has, therefore, two cases to discuss. 

Case 1: If 32 = 0, then, from the first equation in 
Eq. (|5T]) . we have ~ /3'^ ~ 4/i. However, using Eq. ([TT]) 
and Eq. (fT2]) . this result yields a sign contradiction for 
small 1^1 . Thus, this case does not yield any sensible 
solution. 

Case 2: If 51 = 0, then, from the second equation in 
Eq. (|5T]) . we have = (3'^. Therefore, the critical wave 
number for the instability region of the solution Eq. (|14p 
is given by 



(32) 



so that the perturbation functions /i and /2 share the 
same exponential decay rate near ±00 as that of the dark 
soliton uq, where gi and 32 arc constants in the limit of 
I a; I — 00, and the derivatives dgi/dx and d^gi/dx^ are 
bounded in the whole real domain so that 



lim ^ = 0, 

3; I— J- 00 dx 



lim 



(fg^ 
dx^ 



0, 



1,2. 



(29) 



where /? is determined by Eq. ([13]) . It is interesting to 
note that the stability threshold does not depend explic- 
itly on the potential barrier; this is a consequence of em- 
ploying the limit of a; — ±00 in Eq. ([30]) . where the tail of 
the potential decays to zero. Nonetheless, the effect of the 
potential barrier is felt, indirectly, by the DSS through 
the change of its width /3 determined by Eq. ([13]) . 
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